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A  STUDY  OF  WAVENUMBER  INTEGRATION  TECHIQUES 

ABSTRACT 


A  detailed  study  is  made  of  the  Bouchon  (1981)  trapezoidal  integration  rule  for 
evaluation  of  Sommerfeld  integrals.  A  problem  with  non-propagating  arrivab  is 
found  with  integrands  involving  the  zero  order  Bessel  function.  A  mid-point  rec¬ 
tangular  integration  rule  is  offered  as  a  imperfect  way  to  reduce  this  error.  To 
test  numerical  evaluation  of  Hankel  transforms,  the  Haskell  (1963)  wholespace 
solution  is  reformulated,  and  examples  are  given  of  the  analytic,  Bouchon  numeri¬ 
cal  integration  and  mid-point  numerical  integration  of  the  eight  dislocation  and 
two  explosion  Green's  functions. 

Bouchon  (1981)  discussed  the  application  of  a  trapezoidal  numerical  integration  rule  to  the 
evaluation  of  the  Sommerfeld  integral.  The  discussion  followed  previous  work  (Bouchon,  1979)  on 
obtaining  the  solution  for  wave  propagation  due  to  a  point  source  by  evaluation  of  a  two  dimen¬ 
sional  Fourier  transform  over  the  two  spatial  wave  numbers.  Because  the  behavior  of  numerical 
evaluation  of  Fourier  transforms  by  the  discrete  Fourier  transform  technique  is  well  known, 
Bouchon  (1979)  was  able  to  show  that  a  discrete  two-dimensional  trapezoidal  rule  yields  a  wave 
field  corresponding  to  a  distribution  of  point  sources  on  a  rectangular  grid.  Knowing  this,  it  is 
easy  to  establish  the  wave  number  sampling  interval  required  to  yield  seismograms  uncontam¬ 
inated  by  spatial  aliasing. 

Bouchon  (1981)  had  the  objective  of  specifying  the  wave  number  sampling  criteria  to  avoid 
spurious  arrivals  due  to  spatial  and  temporal  aliasing  when  a  trapezoidal  integration  rule  is 
applied  to  the  Sommerfeld  integral.  Bouchon  found  that 

OC 

/F(k)Jo(kr)dk  =  f:<„F(k„)Jo(k.r)Ak  .  (1) 

0  0  =  0 

where  otherwise,  Ak=2»r/L,  ko=nAk,  and  F(k)  is  the  Sommerfeld  kernel.  The 

A 

equality  in  equation  (1)  holds  as  long  as  the  following  two  conditions  hold: 

L  >  2r 

((L-r)*+z^''  >  vt 

where  z  is  the  vertical  distance  between  the  source  and  the  receiver,  r  is  the  radial  distance 
between  the  source  and  receiver,  t  is  the  maximum  time  for  which  a  trace  is  to  be  generated,  and 
V  is  the  velocity  of  the  wave,  of  the  fastest  wave  if  the  problem  has  many  arrivals.  Outside  this 
(r.t)  window  spurious  arrivals  are  seen.  In  addition,  since  a  discrete  Fourier  transform  is  used  to 
invert  (1)  from  the  space-frequency  to  the  space-lime  domain,  Bouchon  used  a  complex  angular 
frequency  given  by  w-ia  to  control  the  inherent  periodicity  in  the  time  series. 

Figures  1-3  illustrate  the  sensitivity  of  the  resultant  seismograms  to  the  parameters  L  and 
Q  The  time  historj'  generated  is  that  of  the  RDS  Green's  function  (APPENDIX),  for  a  point 
source  at  a  depth  of  10  km  beneath  the  receiver.  The  medium  parameters  are  given  in  the  appen¬ 
dix.  Figure  1  has  q=0  0039  and  L=100  km,  Figure  2  has  o=0  03125  and  L  =  100  km,  and  Figure 
3  has  Q=0. 03125  and  L— 200  km.  A  total  of  04  .seconds  of  time  hislorv'  are  presented.  The  effect 
of  the  different  values  of  a  is  not  very  apparent  in  these  figures,  since  the  later  arrivals  are  of  low 
amplitude  to  begin  with.  The  difference  is  seen  in  the  quietness  of  the  traces  prior  to  the  first 
arrival,  at  30  and  70  km  for  example,  and  in  the  occurrence  of  less  ripple  in  Figure  2  than  in  Fig¬ 
ure  1  A  comparison  of  Figures  2  and  3  shows  the  effect  of  increasing  L.  The  number  of  noise 
arrivals  decreases. 

Figure  4  uses  the  sani'-  parameters  as  Figure  3,  except  th.al  a  reduction  velocity  of  6  15 
km 'sec  is  used.  All  traces  start  at  a  lime  t  r  0  15  -  0  50  seconds  .Nl  large  distances,  a 


significant  noise  arrival  overwhelms  the  expected  solution.  This  noise  arrival  appears  to  be  the 
integral  of  the  expected  waveshape  in  the  far-field.  It  is  present  in  Figures  1-3,  and  appears  at  a 
time  corresponding  to  a  non-causal  arrival  traveling  the  vertical  distance  between  the  source  and 
the  receiver.  This  noise  is  seen  only  in  the  integrals  involving  the  Jo(kr)  Bessel  function  and 
corresponds  to  a  k=0  contribution.  The  reason  it  appears  worse  in  Figure  4  than  in  Figure  3,  is 
that  it  wraps  around  to  a  later  time,  when  reduced  travel  times  are  used,  and  is  excessively 
exponentially  increased  when  the  time  series  is  undamped.  This  points  out  the  double  edged  effect 
of  using  complex  frequency,  in  that  the  noise  due  to  later  arrivals  can  always  be  reduced,  but  an 
arrival  earlier  than  the  desired  time  window  will  be  severely  enhanced. 

To  understand  the  problem  and  also  to  appreciate  the  propagating  noise  terms,  we  return  to 
the  Bouchon  (1981)  development.  Bouchon  really  showed  that 

E*.F(k.)J„(k.r)Ak  =  /F(k)J™(kr)dkAk  f)  6(k-kJ  (2a) 

o~0  0  n«s'Oo 

00 

=  /F(k)J„(kr)dk2;r  f)  6(kL-2njr)  (2b) 

0  Oaa-CV> 

CC  ^ 

=  /F(k)J„(kr)dk  S  (2c) 

0  QS-OO 

OO 

=  /F(k)Jn,(kr)dk  (2d) 

0 

OO 

-r  /F(k)J„(kr)dk  5^2cos(nkL) 

0  B-=J 

These  are  essentially  equations  (15-18)  in  Bouchon  (1981),  working  backwards.  We  used  the 
properly  of  the  Dirac  distribution  that  ^(ax)=6(x)/  1  a  j  .  Bouchon  further  expanded  (2d)  to  show 
its  equivalence  to  contributions  of  concentric  rings  of  sources  with  radii  L,  2L,  3L,  ...,  etc,  about 
the  point  source.  We  note  here  that  for  large  kr,  we  can  use  the  asymptotic  expansion  of  the 
Bessel  function  and  a  stationary  phase  approximation  to  show  that  the  left  hand  side  of  (4) 
corresponds  to  an  infinite  set  of  arrivals  in  the  (r,t)  domain  which  are 

sCr.z.il-r  g(L-r,z,t)-i- 1^^  j  g(L+r,z,t)-^...+ 

where  the  function  g  is  the  fhlbert  transform  of  g.  The  validity  of  thb  is  seen  by  examining  the 
noise  arrivals  in  the  Figures  1-4. 

The  functional  form  of  (2d)  is  such  that  the  equation  is  easily  described  in  words!  The  tra¬ 
pezoidal  integration  rule  is  an  approximation  of  the  true  integral,  with  the  second  term  in  (2d) 
being  the  error  term.  .\s  seen  the  error  term  contributes  propagating  numerical  noise.  It  is  also 
apparent  that  the  error  term  is  indeterminate  when  k=0  because  of  the  infinite  summation.  We 
have  numerically  evaluated  the  ten  basic  Green’s  functions  for  dislocation  sources  and  explosive 
sources  (Haskell,  1963;  Haskell,  1964;  Appendix  of  this  paper)  and  found  that  the  k=0  noise 
appeared  only  with  the  integrals  containing  the  Jo(kr)  term.  The  kernel  of  the  Sommerfeld 
integral,  which  is 


-3- 


is  zero  at  k=0,  but  evidently  this  is  not  enough  to  overcome  the  indeterminacy  of  the  error  term. 
On  the  other  hand,  the  negative  radial  derivative  of  (3a). 

does  not  have  such  a  noise  arrival.  The  noise  arrival  is  worst  with  the  negative  vertical  derivative 
of  (3a) 

(3c) 

at  large  distances,  since  the  expected  arrival  decreases  rapidly  due  to  the  radiation  pattern  term 

(— )  while  the  k=0  noise  arrival  does  not. 

K 

A  comparison  of  figures  2  and  3  shows  that  the  k=0  noise  term  is  reduced  in  amplitude  by 
a  factor  of  4,  with  the  same  polarity,  when  Ak  is  decreased  b>  a  factor  of  2,  which  is  expected 
behavior  for  a  trapezoidal  integration  rule  (.Abramowitz  and  Stegun,  1964).  Consider  for  a 
moment  the  trapezoidal  integration  rule,  first  with  sampling  interval  Ak  and  then  with  an  inter¬ 
val  Ak/2.  The  corresponding  rules  are 

CC 

/  f(k)dk  =  Ak(-^f„+f,  +  f2-(-  +  )+0(Ak=^)  (4a) 

0  2 

and 

o 

oc 

/f(k)dk  =  Ak(ifo+f,  ^f,-i-f3+f2+..  +  )+0(Ak*/4)  (4b) 

0  ^22 

where  fn  =  f(nAk).  We  note  that  (4b)  is  just  the  sum  of  (4a)  and  the  mid-point  rectangular 
integration  rule 

00 

/f(k)dk  =  Ak(^  ^fj-t-..  r  )  +  0(Ak-)  (4c) 

0  2  2 

Since  the  numerical  experiment  of  Figures  2  and  3  showed  that  the  k--0  noise  was  reduced  by  a 
factor  of  4  when  Ak  is  reduced  by  a  factor  of  2,  this  suggests  that  the  synthetic  generated  using 

3 

(4c)  must  have  a  k=0  noise  arrival  that  is  —  the  size  of  that  in  (4a)  so  that  (4b)  yields  a  k  =  0 

4 

noise  arrival  that  is  —  that  of  (4a)! 

The  implication  of  thi.s  is  that  a  shifted  rectangular  midpoint  rule  can  be  used  to  substan¬ 
tially  reduce  the  k=0  noise  arrival  In  fact,  we  have  found  through  numerical  experiments  that  a 
numerical  integration  rule 

f:F(kJ.JJk„r)Ak  (5) 

n  -  0 

with 

k„  -  nAk  -  0  218Ak 

works  best  with  Ak  — 27r,  l, 

Figure  5  consists  of  the  same  parameters  as  used  to  generate  Figure  4,  except  that  the 
shifted  rectangular  rule  of  (.'>)  is  used  rather  than  the  trapezoidal  rule  of  (1)  Numerical  noise  is 
still  present  a  low  frequencies,  but  a  least  the  k  0  noise  no  longer  overwhelms  the  expected  sig¬ 
nal  at  large  distance.  The  remaining  propagating  noise  ran  be  reduced  by  using  a  larger  value  for 
L.  .As  indicated  above,  a  phase  change  in  noted  in  the  waveforms  of  the  propagating  noise 
arrivals,  when  they  are  compared  to  the  correspciulmg  arrivals  in  Figure  t 


For  completeness,  it  is  necessary  to  show  that  (5)  is  equivalent  to  the  analytical  integral 
within  the  Bouchon  specified  (r,t)  window.  Following  (2)  we  obtain 

f;F(k.)J„(k.r)Ak  =  /F(k)JJkr)dkAk  £  i{k-kj  (6a) 

B»0  0  ns»-oo 

00 

=  /F(k)J„(kr)dk  (6b) 

0 

~  00 

+  /F(k)J„(kr)dk  5]2cos(n(k-ko)L) 

0  0  —  1 

where  we  define  ka=nAk+ko.  The  only  difference  between  (6b)  and  (2d)  is  the  error  term,  which 
is  no  longer  indeterminate  when  k=0.  For  large  kr,  the  error  term  still  represents  propagating 

arrivals,  with  the  inwardly  and  outwardly  propagating  arrivals  differing  by  phase,  but  now 

the  outwardly  propagating  waves  are  also  phase  shifted  with  respect  to  the  direct  arrivals  from 
the  source. 

The  numerical  integration  rules  used  in  (2)  and  (6)  are  simple  cases  of  general  Newton-Cotes 
integration  rules.  In  numerical  analysis,  one  typically  is  taught  that  a  higher  order  Newton-Cotes 
formula,  such  as  the  Simpson  rule,  yields  better  estimates  of  the  integral.  This  generalization  is 
invalid  when  a  wave  propagation  problem  is  being  solved,  as  is  done  here.  Using  the  same  nota¬ 
tion  as  used  in  (4),  the  Simpson  rule  would  be 

/f(k)dk  =  |^|(f„+4f,+2f2+4f3+...-  )+0(Ak«) 

=  (-j)Ak(-^fo-i-fi+f2i-...+  )  (7) 

+  (y)Ak(fi+f3-t-f5+...-(-  ) 

which  is  recognized  as  a  combination  of  a  trapezoidal  rule  with  sampling  Ak  with  a  midpoint  rec¬ 
tangular  rule  with  sampling  2Ak.  In  evaluating  wave  propagation  integrals  of  the  Sommerfeld 
type,  we  would  see  more  noise  arrivals  that  in  applying  just  the  trapezoidal  rule.  In  thb  case,  an 
attempt  at  additional  numerical  accuracy  backfired.  A  similar  observation  was  made  by  Bakun 
and  Eisenberg  (1970)  in  a  discussion  of  the  numerical  evaluation  of  the  Fourier  transform. 
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APPENDIX:  WHOLESPACE  SOLUTION 

Haskell  (1963)  built  the  solution  for  the  displacement  held  due  to  point  couples  in  a 
wholespace  by  starting  with  the  analytic  solution  for  the  displacement  field  due  to  a  point  force 
given  in  a  cartesian  coordinate  system.  Solutions  for  point  single  couples  were  obtained,  and  the 
solution  was  cast  into  a  cylindrical  coordinate  system,  through  the  use  of  partial  derivatives  of  the 
Sommerfeld  integral.  Haskell  (1964)  extended  the  Haskell  (1963)  work  to  a  layered  halfspace,  to 
include  double-couple  and  dipole  point  forces.  Because  the  Haskell  (1963)  work  gives  both  the 
integrands  of  the  Hankel  transform  as  well  as  the  analytical  answer,  the  wholespace  problem  is 
the  appropriate  one  to  use  for  testing  a  numerical  Hankel  transform  scheme.  The  equations  below 
cast  the  Haskell  (1963)  derivations  into  the  Green’s  functions  for  dislocation  and  explosive 
sources,  given  by  Herrmann  and  Wang  (1985).  The  Green’s  functions  are  defined  as  follow: 


00 


ZDD  =  /F,(k,Lj)Jo(kr)dk 

0 

(la) 

00 

RDD  =  -/F2(k,w)Ji(kr)kdk 

0 

(lb) 

00 

ZDS  =  /F3{k,w)J,(kr)dk 

0 

(Ic) 

00 

RDS  =  /F,(k,u>)Jo(kr)kdk 

0 

(Id) 

00 

-^/;F4(k,u;)  +  F,(k,ui))J,(kr)dk 
f  0 

00 

TDS  =  /F»{k,w)Jo(kr)kdk 

0 

(le) 

00 

4/lF,(k,w)  +  F,(k,w)!J,(kr)dk 
^  0 

00 

ZSS  =  /Fj,(k,ui)J2(kr)dk 

0 

(If) 

00 

RSS  =  /F,(k,w)J,(kr)kdk 

0 

(Ig) 

00 

-^/(F,(k,ui)  ^  F,o(k,w);j2(kr)dk 
f  0 

00 

TSS  =  jFio(k,vj)Ji(kr)kdk 

0 

(Ih) 

00 

/[f'sCk.w)  -f  F,o(k,ui)lJ2(kr)dk 

0 

00 

ZEP  =  /F:(k,w)Jo(kr)dk 

0 

(li) 

00 

REP  =-/F8(k,-;)J,(kr)kdk 

(Ij) 

0 


For  an  arbitrarily  oriented  double  couple  without  moment  source  model  with  vector  n 
=  (nj.no.ns)  normal  to  the  fault  and  f  =  {f1.f2.f3)  if>  the  direction  of  the  dislocation  {  Haskell, 
1963;  Haskell,  1964),  equation  (ll)  of  Wang  and  Herrmann  (1980)  for  the  Fourier  transformed 


displacements  at  the  free  surface  at  a  distance  r  from  the  origin  becomes 

Uj(r,0,'_)  =  ZSS[(fini-f2na)cos2^+(fin2+f2ni)sin2^]  (2a) 

+  ZDS[(fin3+  f3ni)cosiff  +  (f2n3+f3n2)sinv?] 

+  ZDD[f3n3] 

Ur(r.0,!:j)  =  RSS[(fini-f2n2)cos2;?+(fin2+f2ni)sm2^]  (2b) 

+  RDS[(fin3+f3ni)cos^+(f2n3+f3n2)sin'/>] 

+  RDD[f3n3] 

u^(r,0,cj)  =  TSS[{fini-f2n2)sin2y)-(fin2+f2ni)(os2^]  (2c) 

+  TDS[(fin3  +  r3n])sin;p-(f2n3+f3n2)cosi3] 

Explicit  exprc-^sions  for  the  Fj(k,-.')  functions  for  a  point  source  buried  at  a  depth  h  beneath 
the  source  in  a  wliolespacc  with  conipressional  velocity,  a,  shear  velocity,  f),  and  density,  p.  are 
derived  from  Haskell  (1063,  1061)  as  follow: 

Defining 

v^k-  ki  k>k., 
iv/kj^  k* 


\/k=-k/  k>k, 

•>ve  have 


(3a) 


F,(k,w)  =  -^[(2k2-3k2)e'‘'““+3k*e“'n 
4»rpw" 


-"at 


F2(k,w)  ==  — i_[(2k<2-3k2)-^ - +3«/^“'^] 

4jrpw'  t'o 

k2 


FsCk.w)  =  -JL_[2,,„e“'‘^-(2k2-k/)-? - 

inpu"  vg 

F,(k,.;)  =  _±^l2k2e*‘'‘^-(2k2-k^V‘''^] 


Fs(k,w)  = 


4ffpu;' 

k=> 

AtrpJ^ 


k _ ,  k^  -nji 


F.(k,u;)  = 

\i:pu  I'd 


"/Jt, 


F7(k,i*;)  = 


k 


Anpa^ 


Fslk.^)  =  — 

47!fia‘u„ 


(3b) 

(3c) 

(3d) 

(3e) 

(30 

(3g) 

(3h) 

(3i) 

(3j) 


The  vertical  displacement  u,  is  positive  upward,  the  radial  displacement  is  positive  away 
from  the  source,  and  the  tangential  displacement  is  positive  in  a  direction  clockwise  from 

north.  The  vectors  n  and  f  are  still  defined  in  a  local  coordinate  system  at  the  source  in  which 
the  cartesian  axes  are  in  the  north,  east  and  downward  directions.  Following  Herrmann  (1975) 
the  components  of  these  vectors  can  be  expressed  in  terms  of  the  fault  plane  parameters  of  strike, 
dip  and  slip.  The  strike,  <i>i,  is  measured  clockwise  from  north,  the  dip,  df,  is  measured  in  a  posi¬ 
tive  sense  from  the  horizontal  direction  perpendicular  to  strike,  and  the  slip,  Xf,  is  measured  on 
the  fault  plane  in  a  counterclockwise  sense  from  the  horizontal  direction  of  strike.  With  these 
conventions,  all  possible  fault  planes  are  encompassed  by  the  ranges  in  the  angles  of 
0°  <  <?r  <  360°,  0°  <  dr  <  90°,  and  -180°  <  Xf  <  180°.  With  this  notation,  the  sense  of  P-wave 
first  motion  at  the  center  of  the  focal  sphere  is  positive  for  positive  values  of  Xr  and  negative  for 
negative  values.  The  components  of  the  vectors  are 

fi  =  cosXf  cosd>f  -t-  sinXf  cosdf  sin0f 

fo  =  cosXf  sin<if  -  sinX(  cosdf  cosipf 

fs  =  -sinXf  sindf 


nj  =  -sin<^f  sindf 
n;  =  cosi^f  sindf 
n3  =  -cosdf 

Following  Haskell  (1963).  the  analytic  closed  form  solutions  corresponding  to  (la)  to  (Ij)  are 


AnfHjT  dz 


AnfKfT  dz^dr  dt^dv 


AirfHAT  dz^dv  dz^dr 


47rpur  r  drdz  drdz 


47rpar  drdz  dz^ 


dr^  di^dr 


Airpa^  dz 

REP  = 

\vpa‘‘  or 

The  function  Fy  is  the  Sommerfeld  integral, 


-iuiR  CO 

-g-e  =  /— -e“'’^Jo{kr)dk 

n,  Q  l/y 


where 


R*  =  r^  +  h^ 


V*  1.2  W  ‘ 

V  =  k - 


The  specific  expressions  for  the  partial  derivatives  of  the  Sommerfeld  integral 
3Fy  ^fh  ,iW,  h  1 

■ST--'  [f*<T>F 

£El  _  _  L  J  ' 

dr  [R3  ^v^R2 

^Fy  ^ff-9r  15rM  ,iw,f-9r  ISrM  ,iw.of-3r  6r^ 


9h  15h 


9h  15h 


3h  6h 


3h  ISr^z 


1,'VW 
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Following  Wang  and  Herrmann  (1980),  the  following  convention  i  usea  to  define  the 
Fourier  transform  H(f)  of  the  time  series  h{t); 


H(f)  =  j"  h(t)exp(i?rii;  idt 


(5) 


This  integral  is  approximated  by  a  Discrete  Fourier  Transform  in  its  Fast  Fourier  Transform 
implementation  (Brigham,  1974). 

To  generate  synthetic  seismograms,  the  following  source  tina  fii-  ction  is  used: 


275(t)  =  I 


t<0 

0<t<r 

r<t<3T 


--(t/r)-+2(l/r)-l  3r<t<4r 


(6) 


i(t/r)'-4(t/r)+8 


t>4r 


This  time  function  has  a  unit  area.  In  addition,  it  has  spectral  zeros  at  certain  frequencies.  If 
r  =  MAt,  where  M  is  some  power  of  two,  then  spectral  zeros  are  at  frequencies  fjj,  In/S, 

In ^4,  .  .  .  ,  fN/(2M),  where  ffj  is  the  Nyquist  frequency  defined  as  f[^i=  ,  and  N  is  ako  a 

power  of  2.  By  choosing  r  and  At  such  that  one  of  the  spectral  zeros  occurs  at  the  Nyquist  fre¬ 
quency,  the  pulses  can  be  synthesized  and  propagated  through  the  model  without  the  rippling 
introduced  by  an  arbitrar)’,  sharp  high  frequency  spectral  cutoff. 

The  synthetic  seismograms  are  presented  to  show  the  effect  of  using  the  trapezoidal  and 
mid-point  rectangular  numerical  integration  rules.  The  three  sets  of  figures  correspond  to  the 
evaluation  of  (1)  for  a  wholespace.  To  provide  a  uniform  basis  of  comparison,  the  seismic 
moment  is  fixed  at  a  value  of  1.0E-h20  dyne-cm,  the  duration  parameter  r  is  set  to  0.5  seconds, 
and  the  depth  is  fixed  at  a  constant  value  of  10  km.  The  velocities  are  a=6.15km/sec  and 
i=3.55km/sec  and  the  density  is  p=2.8gm/cm*.  A  256  point  time  series  is  synthesized  using 
At=0.25  sec.  Velocity  time  histories  with  units  of  cm/sec  are  generated.  A  causal 
Qo  =  —  10000  is  used,  but  these  are  so  large  that  they  do  not  affect  the  results  displayed.  A 

time  domain  damping  factor  is  used  to  reduce  the  discrete  Fourier  transform  periodicity,  which 
corresponds  to  replacing  all  occurrences  of  w  in  the  frequency  domain  by  w-iO. 046875.  All  resul¬ 
tant  time  series  have  been  undamped. 

.\n  important  aspect  of  the  computations  concerns  the  upper  limit  used  in  the  Hankel 
transform.  Obviously,  a  KNL\X  =  oo  is  out  of  the  question.  Fortunately  the  integrands  become 
small  for  large  k  due  to  the  exponential  decay  terms,  except,  when  the  depth  is  zero.  We  take 


KNIAX  =  max.'F.XCk.  ,-^1 


where  F.XC  is  taken  to  be  3  in  Figures  1-5  and  2  in  the  Figures  in  this  appendix.  The  choice  of 
two  values  controls  special  computations  to  ensure  the  proper  computations  of  the  low  frequency 
contributions  to  the  time  history.  Basically,  an  asym[>totic  expansion  is  used  in  this  case  to  deter¬ 
mine  the  contribution  from  k=KNL‘XX  to  k— oo.  Too  large  a  value  for  FAC  will  require  too 
many  computations  without  much  discernible  difference  in  the  high  frequency  synthetics,  while 
too  low  a  value  will  introduce  significant  low  frequency  errors  in  the  time  series. 


NUMERICAL  EXAMPLES 

The  first  set  of  figures,  "whom”  gives  the  analytic  solution  using  equation  (4).  A  reduced 
travel  time  plot  with  initial  time  given  by  t  =  -1.01  -i-  r/6.15  seconds,  where  r  is  the  epicentral 
distance.  The  correspondence  between  the  identifier  JSRC  and  the  specific  Green’s  function  is  as 
follows: 

JSRC  Green’s  Function 

1  ZDD 

2  ROD 

3  ZDS 

4  RDS 

5  TDS 

6  ZSS 

7  RSS 

8  TSS 

9  ZEP 

10  REP 

The  Green’s  functions  show  the  P-wave  and  S-wave  arrivals  expected  at  large  distances.  The  S- 
wave  arrival  on  the  RDD  and  RSS  components  has  a  waveform  that  is  an  integral  of  the  expected 
far-field  arrival,  the  corresponding  shape  of  the  P-wave. 

The  second  set  of  figures,  "WlOOOTO. 00000”,  is  the  result  of  the  numerical  evaluation  of  the 
Hankel  transforms  using  the  Bouchon  trapezoidal  rule.  An  L  =  1000  km  was  used.  The  k=0  noise 
is  very  apparent  on  the  ZDD,  RDS,  TDS  and  ZEP  traces.  This  noise  arrival  also  demonstrates  the 
periodicity  of  the  discrete  Fourier  transform  as  well  as  the  problems  with  using  complex  fre¬ 
quency.  Propagating  noise  arrivals  are  seen  in  the  traces  beyond  300  km,  because  the  criteria 
relating  L,  r,  z,  v  and  t^j^,  in  equation  (1)  are  no  longer  satisfied.  To  eliminate  these,  we  need  only- 
make  L  somewhat  larger.  Nothing  of  value  should  be  expected  of  the  trace  at  500  km  because 
the  direct  and  first  inwardly  propagating  noise  arrival  superimpose.  At  this  distance,  the  inequali¬ 
ties  of  the  Bouchon  analysis,  equation  (l),  are  violated. 

the  last  set  of  figures,  "WlOOOTO.21739”  ,  corresponds  to  the  use  of  the  mid-point  rule  (5) 
with  ko  =  0.21739Zik.  As  designed,  the  k=0  noise  arrival  is  significantly  reduced,  although  some 
low  frequency  numerical  noise  is  introduced  at  large  distance.  The  difference  in  the  integration 
rules  is  most  readily  apparent  in  the  TDS  and  ZEP  time  histories. 


FIGURE  CAPTIONS 

Fig.  1.  Synthetic  sebmograms  for  the  RDS  component,  using  o=0.0039  and  L=100km. 

Fig.  2.  Synthetic  sebmograms  for  the  RDS  component,  using  q=0.03125  and  L=  100km. 

Fig.  3.  Synthetic  sebmograms  for  the  RDS  component,  using  o=0.03125  and  L=200km. 

Fig.  4.  Synthetic  sebmograms  for  the  RDS  component,  using  a=0.03125  and  L=200kn,  but 
using  a  reduced  travel  ti:»ie  display. 

Fig.  5.  Synthetic  seismograms  for  the  RDS  component,  using  a=0.03125  and  L=200km,  a 
reduced  travel  time  display,  but  using  the  shifted  rectangular  integration  rule  rather 
than  the  Bouchon  trapezoidal  integration  rule. 
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3.  Synthetic  seismograms  for  the  RDS  component,  using  q=0. 03125  and  L=200km. 
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Fig  5.  Synthetic  seismograms  for  the  HDS  component,  using  q— 0.03125  and  L  -200km,  a 
reduced  travel  time  display,  but  using  the  shifted  rectangular  integration  rule  rather 
than  the  Bouchon  trapezoidal  integration  rule. 

JSRC=H  RELATIVE  AMPLITUDE 


The  following  figures  provide  the  ten  Green’s  functions  for 
the  analytic  solution,  "whom,"  listed  at  the  bottom  of  the  model 
page,  for  the  Bouchon  integration  scheme  "W lOOOTO. 00000 ," 
and  for  the  modified  rectangular  rule,  "W1000T0.21739 

Note  that  the  modified  rectangular  rule  does  reduce  the  false 
k  =  0  arrival.  However,  the  synthetics  obtained  using  the  numerical 
integration  techniques  still  have  noise  arrivals  at  large  distances 
due  to  too  small  a  choice  for  L.  This  noise  is  seen  at  distances 
larger  thein  300  km  and  corresponds  to  a  P-wave  arrival  at  a  time 
of  (1000  -  r)/6.15  seconds.  This  is  readily  seen  in  the  JSRC  =  2,  9 
or  10  Green's  functions.  There  is  no  such  problem  with  S  arrivals, 
as  seen  in  the  predominantly  SH  Green’s  functions  JSRC  =  5  and  8. 
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JST^C=9  RELATIVE  AMPLITUDE 
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■1.  108E-06 

100.0000 


959E-07 

150.0000 


2. 592E-07 

200. 0000 


4.785E-07 
250. 0000 


-3. 828E-06 
500. 0000 


AMPLITUDE 


37 


-1.117E-05 

1.0000 


-H.  H6HE-06 
50. 0000 


V 


-2.  109E-06 

100. 0000 


-1.HH2E-06 

150.0000 


-1.087E-06 

200.0000 


-2.693E-07 
250. 0000 


-7.  255E-07 
300. 0000 


-6.188E-07 
350. 0000 


JSPC=7  RELATIVE  AMPLITUDE 


ALPHA =04 0H7 
DEPTH=10. 000 
FL  =0*000 
FU  =2.000 
DT  =0.250 

Nr  Nl.  N2l256rlrl29 


D 

A 

B 

HHO 

QA  INV 

11.000 

6.150 

3.550 

2.800 

0.00010 

6.150 

3.550 

2.800 

0.00010 

ND  =11 

DMIN=1.00000E+02 
DD  =0. OOOOOE+00 


TMIN  =0.0000 
TMAX  =63.75000 
CRrTD 


QB  INV 

0.00010 

0.00010 


PABABOLIC  ITYPE=H 
WIOOOTO. 21739 


2 


3.363E-0S 

KOOOO 


50.0000 


2.  066E-06 


-l.*l29E-06 
ISO. 0000 


-1.082E-06 

200.0000 


<e.662E-07 

250.0000 


-7.  229E-07 
300.0000 


^ 196E-07 
350.0000 


5.  ‘»31E-07 
H00.0000 


H.  853E-07 
^50.0000 


1.596E-07 
500.  0000 


JST^C=2  l^ELATIVE  AMPLITUDE 


1.393E-05 

KOOOO 


5* 12HE-07 
150«0000 


^.6H6E-07 

200.0000 


-1.715E-07 

250.0000 


JSRC=H  RELATIVE  AMPLITUDE 


6 


4.125E-06 

1*0000 


■^.500E-06 

50*0000 


■l*118E-06 
100* 0000 


5*0‘t2E-07 
150* 0000 


2*599E-07 
200*  0000 


■l*705E-07 

250*0000 


•l*2H6E-07 
300* 0000 


9*2H3E-08 

350*0000 


7*303E-08 
HOO* 0000 


6*3H6E-08 
H50* 0000 


1* 106E-07 
500* 0000 


JSRC=6  RELATIVE  AMPLITUDE 


9 


•1.  128E>0S 

1.0000 


■2.251E-05 
SO. 0000 


■1.  108E-05 

100. 0000 


7.  *»38E-06 
150. 0000 


5. 107E-06 

200. 0000 


4. 25HE-06 
250. 0000 


'3.701E-06 
300. 0000 


'3.  158E-06 
350. 0000 


2. 6H5E-06 
HOO. 0000 


TIELATIVE  AMPLITUDE 


?  /  i 


51 


2.52SE-06 

UOOOO 


^•061E-04 

50*0000 


0H8E-06 

100*0000 


•1*  HOSE -06 
150*0000 


•l«06SE-06 

200*0000 


■e*5H8E-07 

250*0000 


7*2H8E-07 

300*0000 


6«112E-07 
350*  0000 


-5*3H8E-07 

H00«0000 


■H*753E-07 

H50*0000 


.  I TUBE 


